Carga paquetes
library(sf)
library(spdep)
library(tmap)
library(INLA)
# library(mapview)
# library(leaflet)
# library(leafem)
# library(sp)
#
#
# library(ggplot2)
# library(ggpubr)
# library(knitr)Al comenzar la pandemia producida por COVID-19, surgió el interés de analizar la distribución espacial de los riesgos infección.
Ajustamos un modelos lineal generalizado con efectos aleatorios espaciales con INLA para estimar riesgo de aparición del virus en diferentes radios censales de la Ciudad de Córdoba, Argentina.
Se trabajará con una base de datos vectorial (.gpkg) que incluye polígonos de los radios censales de Córdoba y Gran Córdoba. Para cada radio censal se cuenta con información de
# Generación de lista con vecindarios
nb <- poly2nb(datos)
head(nb)
plot(st_geometry(datos), border = "grey", lwd = 0.5)
plot(
nb,
coordinates(as(datos, "Spatial")),
add = TRUE,
col = "blue",
points = FALSE,
lwd = 0.5
)
## [[1]]
## [1] 2 3 4 29 37 237 251 252 258 266 267 1498 1500
##
## [[2]]
## [1] 1 3 8 14 37
##
## [[3]]
## [1] 1 2 4 8
##
## [[4]]
## [1] 1 3 5 6 7 8 266 267
##
## [[5]]
## [1] 4 6 267
##
## [[6]]
## [1] 4 5 7 10 11 270
# Ajuste del Modelo inflado en ceros
summary(datos)
formula <-
Casos ~
Fragmentacion +
HogaresNBI +
HogaresHacinamiento +
HogaresJefe.Univ. +
Bancos +
f(re_u, model = "besag", graph = g) +
f(re_v, model = "iid")
# formula <-
# Casos ~ Bancos + ValorTierra + f(re_u, model = "besag", graph = g) + f(re_v, model = "iid")
res_inflpoi <-
inla(
formula,
family = "zeroinflatedpoisson1",
data = datos,
E = E,
control.predictor = list(compute = TRUE)
)
summary(res_inflpoi)## Poblacion Fragmentacion ValorTierra HogaresNBI
## Min. : 7.0 Min. :1.000 Min. : 0.01171 Min. : 0.0000
## 1st Qu.: 629.8 1st Qu.:2.000 1st Qu.: 0.14349 1st Qu.: 0.3755
## Median : 835.0 Median :4.000 Median : 0.35141 Median : 1.0855
## Mean : 873.2 Mean :3.367 Mean : 1.00000 Mean : 1.6939
## 3rd Qu.:1086.0 3rd Qu.:4.000 3rd Qu.: 0.73210 3rd Qu.: 2.4040
## Max. :3273.0 Max. :4.000 Max. :12.88503 Max. :20.6897
## HogaresHacinamiento HogaresJefe.Univ. HogaresConComput. HabSup
## Min. :0.0000 Min. : 0.0000 Min. : 0.00 Min. : 0.006
## 1st Qu.:0.1364 1st Qu.: 0.9802 1st Qu.:11.96 1st Qu.: 44.619
## Median :0.5051 Median : 3.3943 Median :18.05 Median : 73.769
## Mean :0.8617 Mean : 5.5653 Mean :20.02 Mean : 83.815
## 3rd Qu.:1.2639 3rd Qu.: 9.0126 3rd Qu.:25.26 3rd Qu.: 97.168
## Max. :6.6606 Max. :26.6055 Max. :57.17 Max. :606.024
## Bancos Farmacias EstSaludInt Ban_Farm_Salud
## Min. :0.00000 Min. :0.0000 Min. :0.0000 Min. :0.0000
## 1st Qu.:0.00000 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:0.0000
## Median :0.00000 Median :0.0000 Median :0.0000 Median :0.0000
## Mean :0.05301 Mean :0.2446 Mean :0.1536 Mean :0.4512
## 3rd Qu.:0.00000 3rd Qu.:0.0000 3rd Qu.:0.0000 3rd Qu.:1.0000
## Max. :8.00000 Max. :6.0000 Max. :6.0000 Max. :9.0000
## Casos geom E SIR
## Min. : 0.0000 MULTIPOLYGON :1660 Min. :0.002241 Min. : 0.000
## 1st Qu.: 0.0000 epsg:22174 : 0 1st Qu.:0.201595 1st Qu.: 0.000
## Median : 0.0000 +proj=tmer...: 0 Median :0.267299 Median : 0.000
## Mean : 0.2795 Mean :0.279518 Mean : 1.007
## 3rd Qu.: 0.0000 3rd Qu.:0.347649 3rd Qu.: 0.000
## Max. :36.0000 Max. :1.047749 Max. :89.967
## re_u re_v
## Min. : 1.0 Min. : 1.0
## 1st Qu.: 415.8 1st Qu.: 415.8
## Median : 830.5 Median : 830.5
## Mean : 830.5 Mean : 830.5
## 3rd Qu.:1245.2 3rd Qu.:1245.2
## Max. :1660.0 Max. :1660.0
##
## Call:
## c("inla.core(formula = formula, family = family, contrasts = contrasts,
## ", " data = data, quantiles = quantiles, E = E, offset = offset, ", "
## scale = scale, weights = weights, Ntrials = Ntrials, strata = strata,
## ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose =
## verbose, ", " lincomb = lincomb, selection = selection, control.compute
## = control.compute, ", " control.predictor = control.predictor,
## control.family = control.family, ", " control.inla = control.inla,
## control.fixed = control.fixed, ", " control.mode = control.mode,
## control.expert = control.expert, ", " control.hazard = control.hazard,
## control.lincomb = control.lincomb, ", " control.update =
## control.update, control.lp.scale = control.lp.scale, ", "
## control.pardiso = control.pardiso, only.hyperparam = only.hyperparam,
## ", " inla.call = inla.call, inla.arg = inla.arg, num.threads =
## num.threads, ", " blas.num.threads = blas.num.threads, keep = keep,
## working.directory = working.directory, ", " silent = silent, inla.mode
## = inla.mode, safe = FALSE, debug = debug, ", " .parent.frame =
## .parent.frame)")
## Time used:
## Pre = 2.28, Running = 36.1, Post = 1.44, Total = 39.8
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) -0.662 0.492 -1.641 -0.658 0.287 -0.648 0
## Fragmentacion 0.169 0.106 -0.036 0.168 0.378 0.167 0
## HogaresNBI -0.052 0.070 -0.199 -0.050 0.078 -0.044 0
## HogaresHacinamiento 0.115 0.140 -0.154 0.113 0.397 0.108 0
## HogaresJefe.Univ. 0.052 0.023 0.008 0.052 0.097 0.052 0
## Bancos 0.138 0.194 -0.249 0.140 0.514 0.145 0
##
## Random effects:
## Name Model
## re_u Besags ICAR model
## re_v IID model
##
## Model hyperparameters:
## mean sd 0.025quant
## zero-probability parameter for zero-inflated poisson_1 0.598 0.052 0.503
## Precision for re_u 0.346 0.062 0.242
## Precision for re_v 31.850 40.820 5.001
## 0.5quant 0.975quant
## zero-probability parameter for zero-inflated poisson_1 0.595 0.706
## Precision for re_u 0.340 0.485
## Precision for re_v 19.933 134.836
## mode
## zero-probability parameter for zero-inflated poisson_1 0.581
## Precision for re_u 0.327
## Precision for re_v 10.215
##
## Marginal log-Likelihood: -2315.06
## is computed
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
datos$RR <- res_inflpoi$summary.fitted.values[, "mean"]
datos$RR_LI <- res_inflpoi$summary.fitted.values[, "0.025quant"]
datos$RR_LS <- res_inflpoi$summary.fitted.values[, "0.975quant"]
popup_vars <- c(
"Fragmentación (Nivel)" = 'Fragmentacion',
"Valor Tierra Medio ($/m2)" = 'ValorTierra',
"Hogares NBI (%)" = 'HogaresNBI',
"Hogares Hacinamiento (%)" = 'HogaresHacinamiento',
"Hogares Jefe Univ.(%)" = 'HogaresJefe.Univ.',
"Bancos",
"Casos",
"Población" = 'Poblacion',
"E",
"SIR",
"RR_LI",
"RR",
"RR_LS"
)tmap_mode('view')
mapas <-
tm_basemap(c(Urbano = "OpenStreetMap", Satelite = "Esri.WorldImagery")) +
tm_shape(datos,
name = 'Casos') +
tm_polygons(
col = "Casos",
id = "Casos",
border.col = "gray50",
border.alpha = .5,
style = "fixed",
title = "Casos por Radio",
palette = "YlOrBr",
breaks = c(0, 1, 2, 3, 5, 10, 15, 25, 36),
popup.vars = popup_vars,
legend.format = list(
scientific = TRUE,
format = "f",
digits = 0
)
) +
tm_shape(datos,
name = 'Tasa de Infección - SIR') +
tm_polygons(
col = "SIR",
id = "SIR",
border.col = "gray50",
border.alpha = .5,
style = "fixed",
palette = hcl.colors(7, "ag_GrnYl"),
breaks = c(0, 5, 10, 15, 20, 25, 30, 50, 60, 90),
popup.vars = popup_vars,
legend.format = list(
scientific = TRUE,
format = "f",
digits = 0
)
) +
tm_shape(datos,
name = 'Riesgo Relativo') +
tm_polygons(
col = "RR",
id = "RR",
style = "fixed",
palette = "viridis",
alpha = 0.8,
title = "Riesgo Relativo",
breaks = c(0.4, 1, 2, 4, 8, 10, 15, 20, 30, 50, 85, 96),
popup.vars = popup_vars,
border.col = "gray50",
border.alpha = .5,
legend.format = list(
scientific = TRUE,
format = "f",
digits = 1
)
)
mapas